Simple strong glass forming models: mean— field solution with activation 



Arnaud Buhot 1 ' 2 *, Juan P. Garrahan 1 ^ and David Sherrington 1 * 
1 Theoretical Physics, University of Oxford, 1 Keble Road, Oxford, 0X1 SNP, U.K. 
2 CEA Grenoble, DRFMC, SI3M, 17 rue des Martyrs, 38054 Grenoble, France. 
C<\ ' (February 1, 2008) 

o : 

We introduce simple models, inspired by previous models for froths and covalent glasses, with 
. trivial equilibrium properties but dynamical behaviour characteristic of strong glass forming systems. 

These models are also a generalization of backgammon or urn models to a non-constant number 
of particles, where entropic barriers are replaced by energy barriers, allowing for the existence of 
' activated processes. We formulate a mean-field version of the models, which keeps most of the 

features of the finite dimensional ones, and solve analytically the out-of-equilibrium dynamics in 
the low temperature regime where activation plays an essential role. 
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I. INTRODUCTION 



Glass forming systems, like supercooled liquids, display remarkable dynamical properties, and are the subject of 
intense experimental and theoretical studies (for reviews see [1-3]). Out of all glass formers, the ones which have 
the seemingly simpler macroscopic behaviour are strong liquids, typically covalently bonded systems like Si02 and 
Ge02, in which viscosities and relaxation times increase exponentially with decreasing temperature, following a simple 
Arrhcnius law [1]. While the increase of timescales in strong systems is not as dramatic as in fragile ones, it is still 
i-^j j very pronounced, and, moreover, it is not accompanied by the growth of static lengthscales or similar corresponding 
c"j . structural signatures. 

The separation of statics and dynamics is a central feature of structural glass formers [4] , and is in marked contrast 
i with other dynamically arrested systems, like spin glasses and other systems with quenched disorder [5,6]. This leads 

naturally to an approach for modeling these systems based on the idea that glassiness is not a consequence of disorder 
t-H , or frustration in the static interactions but of constraints on their dynamics [7-10]. The aim of this paper is to build 
on our earlier efforts in the modeling of covalent systems and froths [11—13], where the dynamical constraints that lead 
to glassy slowdown can be related directly to elementary structural processes, and develop a class of models which 
' retain the essential features but are simple enough to allow for analytical solutions, and study in detail the connection 
with underlying effective diffusion-annihilation processes [14,15] which we believe are an essential component of strong 
' glass formers in general. 

The models we study in this work are described explicitly and concisely in section II. Here we outline briefly the 
background considerations which were their conceptual basis. They are a distillation of a series of minimalist models 
exhibiting strong glass-like dynamical behaviour through the operation of kinetic constraints despite non-interacting 
Hamiltonians and trivial thermodynamics (with no equilibrium phase transitions) and no imposed quenched disorder. 

The first model [11] was a 2-dimensional topological 'foam' consisting of a fully pairwise connected network of 
three-armed vertices with energy function E — — 6) 2 where 7ij is the number of edges (sides) of cell i ; thus it 

emulates desire for local hexagonal structure without the complication of interaction. The dynamics is Tl (see Fig. 1) 
performed stochastically with acceptance probabilities determined by the Boltzmann ratio exp(— AE/T) where AE is 
the resultant energy change and T is the temperature. This ensures eventual equilibration at any finite temperature 
but also provides a kinetic constraint since each Tl move involves four cells, decreases n by 1 for each of the two 
initially adjacent cells and increases n of the other two by the same amount. 
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FIG. 1. A Tl move. The initially adjacent cells decrease their number of neighbours by 1 whereas the cells which become 
adjacent increase their number of neighbours by 1. 



This model can be considered variously as an idealization of a glass of covalent sp 2 hybrid bonds which change 
connections, or of the Voronoi cells surrounding a fluid of atoms, or of a foam ; in the first case effectively allowing for 
harmonic changes of vertex angle energy but ignoring cell-cell correlation energies, in the second and third replacing 
continuous energies associated with interatomic forces or surface tension by the simple topological energies. At low 
temperatures the model exhibits two-timescale macro-behaviour, for example in the decay of the energy from a 
random start or in the autocorrelation function even in equilibrium. The fast timescale is temperature independent 
and the slow timescale is Arrhenius-like: r ~ exp(a/T), a ~ 2 [11,12]. This behaviour can be understood by 
considering m = 6 as 'perfect' and m ^ 6 as 'defect' cells, and the cooperative behaviour in terms of the diffusion and 
annihilation of defects. The fast processes correspond to ones which involve either an energy decrease or no energy 
change, while the slow processes involve energetic defect creation. Defects with \m — 6| > 1 disappear rapidly at low 
temperature so we concentrate on \m — 6| < 1. Let us denote n = 5(7) cells by A(A) and n — 6 by 0. Then energy 
reduction is due to the annihilation of an AA neighbouring pair (or dimer). The Tl process however involves four 
cells and only allows A A dimer annihilation if at most only one of the other cells involved is a 0; thus the annihilation 
of two dimers [16] 

2A + 2A^40 (1) 

exists as well as the single dimer annihilation [17] 

2A + A + 30 + A, 2A + A + -> 30 + A. (2) 

AA dimers can diffuse on the fast (microscopic attempt) timescale if the other pair of cells in the Tl process are both 
0: 

A + A + 20^ 2Q + A + A. (3) 

Defects are only annihilated in dimer pairs (AA) and isolated defects must be paired to annihilate, hence they must 
be moved together. However they can only move via processes which involve the creation of an AA dimer and a 
consequent energetic penalty; i.e. by the inverses of processes (2) 

A + 30 -» + 2A + A, A + 30 -» + 2A + A, (4) 

which are slower by the Boltzmann factor exp(— 2/T). Dimer pairs, of either of the two possible AA combinations 
on the right hand sides of processes (4) can then diffuse away freely by process (3). Processes (2) themselves also 
allow slow movement of A (or A) via collision with a roaming AA dimer, leading to annihilation of all three A, A 
but creation of A (or A) on the cell which was originally 0. This also involves a slowing factor exp(— 2/T) since the 
equilibrium density of AA dimers involves such a factor. As a consequence, at low temperatures, the energy after a 
quench from a random/high temperature start quickly decays to a plateau through elimination of AA dimers, leaving 
mainly (and at T = only) isolated defects and clusters of only like defects. After this plateau, the energy further 
decays to equilibrium with the Arrhenius-like timescale r ~ exp(a/T) due to the effective dimer-mediated diffusion 
of isolated defects and eventual AA pairing and annihilation. 

The second generation of models is built on the above picture, replacing the foam by lattice-based models with 
'spins' associated with the dual plaqucttes [18] 

Si = 1,0,-1, (5) 

the energy by 

T=jf>f, (6) 

i=l 
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and the Tl move by one involving four spins around a corresponding 'Feyman diagram' edge in a hexagonal lattice, 
a square vertex in a square lattice, etc. [13]. For J > this behaves essentially as for the foam model [19], but allows 
much larger simulations. The resultant better-accuracy data permits comparison with the asymptotic predictions of 
the field theory of simple diffusion-annihilation processes [14,15] of type A + B — > (in the usual notation of that 
subject); the agreement is quite good for both decays which can be viewed as diffusion-annihilation processes with 
different timescales [13]. 

The second class of models permits a further modification, the taking of J < 0. In this case the ground state is 
highly degenerate (Si = ±1 independently on each plaquette), as compared with the unique ground state for J > 
(Si = 0; all i), but the defects are of only one type (Si = 0). The two timescale behaviour of the energy decay to 
equilibrium and of the equilibrium autocorrelation functions persists. Ignoring the complications of the ground state 
degeneracy and further constraints imposed by the Tl process (\Si\ cannot be increased above 1) the behaviour is 
now describable in terms of only A, dynamics (the now referring to S = ±1 and A to S = 0). The qualitative 
discussion earlier still applies with A (or B) replaced by A. Again the comparison with the predictions of simple field 
theory of diffusion-annihilation is fair, but less precise, presumably in part due to the more complicated ground state 
degeneracy and constraints. 

The models studied in this paper are a further simplification which maintain the features of (i) fast annihilation 
of dimers, (ii) fast diffusion of dimers, (iii) motion of isolated defects only by slow creation of dimers. It allows for 
cither a single defect type (A) and identical defect dimers (AA) or two defect types (A and B) and mixed dimers 
(AB), but in both cases with a non-degenerate (T = absorbing) ground state (0). They also allow separation of 
the timescales for processes (i) and (ii) above, as well as their separation from (iii). The explicit formulation is given 
in the next section, but the above should clarify why they are appropriate. 

The rest of the paper is organized as follows. In section II we give a full description of the models studied in this 
work. Their general features are discussed and compared with those of previous models is section III. In section IV 
we solve analytically the mean-field version of the models. Our conclusions are given in section V. Details of the 
analytical calculations are provided in the appendices. 



II. DESCRIPTION OF THE MODELS 



The minimalist models we study here are based on a coarse-grained simplification of the ideas described above. 
They correspond to defects (or particles), which can be either of a single kind A, or different kinds, A and B, and 
live in a d dimensional lattice. They can also be considered as a generalization of backgammon [21] or urn models 
[22] to a non-constant number of particles, with energetic barriers rather than entropic ones, thus allowing for the 
existence of activated processes. The analogy with previous models [11-13] is obtained through the dynamical rules 
which mimic the different processes (2), (3) and (4). The two models under consideration are the following. 



A. Single type of particles 

As for previous models, with each defect (particle in our case) we associate a unit energy ( J = 1) leading to the 
Hamiltonian 

N 

ff = £«i, (7) 

i=l 

with Hi the occupation number on site i. Due to the non- interacting Hamiltonian, the equilibrium properties are 
trivial. The probabilities to have n particles on a given site at temperature T — 1/0 are: 

P^(0)=e-^[Y,e-A (8) 

\n=0 / 

and lead to an equilibrium energy or concentration of particles: 

c eq (0) = £ np%(0) (9) 

n=0 
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with n max the maximum number of particles per site. In the following we use the smallest number (n max = 3) 
compatible with the dynamical rules. The infinite temperature concentration c eq (j3 = 0) = 3/2 corresponds to equal 
probabilities = 1/4 whereas the low temperature concentration vanishes as c eq ((3) ~ e~^. 

The dynamical rules are inspired from the Tl moves [11-13]. Three different kinds of moves are considered [20]: 

(i) The annihilation of two particles analogous to processes (2): three particles disappear from site i and only one 
appears on a neighbouring site j with a rate 1 

(rii, rij) — > (rii - 3, rij + 1) (10) 

(ii) The dimcr diffusion analogous to process (3): two particles move from site j to a neighbouring site j with a 
diffusive rate D 

rij) — > (rii - 2, rij + 2) (11) 

(iii) The creation of two particles analogous to processes (4): a particle disappears from site i to create three particles 
on a neighboring site j with a rate e~ 2/3 

(rii, rij) — > (rii - 1, + 3) (12) 

Those moves have also to respect the maximal number of particles per site n max = 3. The rates considered satisfy 
detailed balance with respect to Hamiltonian (7), ensuring that the equilibrium properties [equations (8) and (9)] 
will be reached asymptotically by the dynamics. All these processes are of the form (rii, n j) — * (n% — x, rij + y) with 
x + y = 4, which reflects the four cell character of the original model transitions. 

The introduction of a diffusive constant D allows to separate explicitly the timescales for diffusion of dimers from 
that of the annihilation process. It is also interesting to notice that those rules are defined for any dimension and any 
lattice, even for a fully connected network, leading the possibility to study the mean-field version of the model. 



B. Two different types of particles 

A possible modification of the model consists in considering two different types of particles (A and B) to keep the 
analogy with the possible sign of the topological charges qi = 6 — rii of the cells in the foam model. As a consequence, 
the dynamical rules are closer to the Tl move. This modification allows to stress the relation between the last decay 
of the energy or concentration of particles and diffusion-annihilation processes [14,15], since the critical dimension 
and exponents for the A + B — > process is different from those of the A + A — > one. 

The presence of two types of particles introduce modifications to the equilibrium properties. In the following we still 
consider a restriction on the number of particles per site (n max = 3) irrespective of their nature, but, in addition, the 
difference between the numbers of A and B particles on a site is limited to —1, or 1 (this last restriction absent within 
the previous model avoids configurations with all A's or all B's on a given site). As a consequence the equilibrium 
probabilities to have ha particles A and ns particles fiona given site are (with n = riA + ris)'- 

PZ,n B W - |e-^6(n max - n)Q(l - \n A - n B \) (13) 
leading to the equilibrium concentration of particles: 

c eq (f3)= ]T npZ, nB (j3) (14) 

TlA,nB 

with: 

Z = E ^"©(^max-^ea-K-nsI) (15) 
n A ,n B 

and 0(x) the Heaviside function (Q(x) — 1 if x > and Q(x) — otherwise). 

The infinite temperature concentration of particles is c eq ((3 = 0) = 5/3 [six possible configurations with equal 
probabilities: (riA,riB) — (0, 0), (1,0), (0, 1), (1, 1), (2, 1) and (1, 2)] whereas for low temperatures the concentration 
vanishes as c eq (/3) ~ 2e~' 3 (the coefficient 2 is coming from an entropic effect due to the two types of particles). 
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The dynamical rules are a straightforward generalization of the Tl moves [11-13]. Three different kinds of moves 
are considered: 

(i) annihilation of an AB dimer analogous to processes (2), an AB dimer and another particle disappear from site i 
but only the single particle appears on a neighbouring site j with a rate 1 

[(AAB)i, (X),} [(0),, (AX)j], [(ABB),, (X),} [(0),, (BX)A (17) 

(ii) AB dimer diffusion analogous to process (3), a dimer moves from site i to a neighbouring site j with a diffusive 
rate D 

[(ABX) l ,(Y) 3 ]^[(X) t ,(ABY) ] ] (18) 

(iii) creation of an AB dimer analogous to processes (4), a single particle from site i move to a neighbouring site j 
creating a dimer on this site with a rate e~ 2/3 

[(AX) U (0),] [(X) u (AAB)j], \(BX) U (0),] [(X) u (ABB),} (20) 

In all of these processes, symbols X and Y stand for possible A, B or particles respecting the restrictions in the 
number of particles on each site. The rates again satisfy detailed balance conditions, ensuring equilibration. 

III. GENERAL FEATURES 

We now show that the two models introduced in the last section share common behaviour for all dimensions and with 
the models considered previously [11-13]. In what follows we discuss equilibrium dynamical properties, in particular 
the existence of two different timescales, as well as out-of-equilibrium features like the multi-stage decay of the energy 
density after a quench. 

A. Dynamics in equilibrium 

Let us first consider the equilibrium dynamical properties of these systems. As we already mentioned, the dynamical 
behaviour is characteristic of a strong glass [1]. The relaxation time increases exponentially with the inverse tempera- 
ture, corresponding to an Arrhenius law. A simple way to determine this relaxation time is from the auto-correlation 
function: 

C(t,t') = (mMmV)) (21) 

with the brackets denoting ensemble average. In equilibrium this two-time function reduces to a single time equilibrium 
correlation C eq (t — t') due to the time translational invariance. We can define the relaxation time r from C^ q (r) = 
Cg (0)/e, where the connected correlation C° (t) — C eq (t) — c 2 eq . At low temperatures and for all diffusive constants, 
the temperature dependent and slower process is the creation of particles which has energy barrier AE = 2. As a 
consequence we expect the following Arrhenius law for the relaxation time: 

r(/3)oce 2/3 . (22) 

This behaviour is confirmed by numerical simulations. The expected Arrhenius law (22) is recovered for all dimensions 
and all diffusive constants (see Fig. 2). Similar results are found for the model with two types of particles. (Simulations 
have been performed using continuous time Monte Carlo [23] for systems of N = 10 6 sites.) 
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Inverse Temperature 

FIG. 2. Relaxation time as function of the inverse temperature for different dimensions (d — 1,2,4 and oo) and a diffusive 
constant D = 10~ 4 . The line corresponds to the expected r oc e 2/3 behaviour. 

The equilibrium correlation explicitly shows the two timescale behaviour (see Fig. 3): (i) dimer diffusion or annihi- 
lation on a short timescale independent of the temperature and (ii) isolated particle motion through the creation of 
a dimer on the relaxation timescale. At low temperatures when the two timescales are well separated the correlation 
presents a plateau separating the two relaxing decays. The relative position of this plateau is increasing with decreas- 
ing the temperature due to the fast disappearance of dimers. The decay to zero of the connected correlation is also 
dimension dependent. A change of behaviour occurs between the mean-field case (d — oo), where it is exponential, 
and finite dimensions (d = 2 in Fig. 3), where it only vanishes algebraically. The probability for a particle to come 
back to the same place (which depends on the dimension of the system) explains this difference. 




Time 

FIG. 3. Normalized equilibrium autocorrelation Cl q {t) / Cl q {Q) for the model with a single type of particles and a diffusive 
constant D = 1. Different temperatures are considered (from left to right f3 — 1, 2, 3, 4, 5 and 6) as well as dimensions: d = 2 
(dashed curves) and mean-field d — oo (full curves). 



B. Out— of— equilibrium dynamics 

We now consider the out-of-equilibrium behaviour of the models, in particular the decay of the concentration of 
particles (energy density) c(t) = N (H (t)) , after a quench from an infinite temperature to a low temperature T at 
time t = 0. This decay shows an interesting structure with two intermediate plateaux when the diffusive constant D 
and the final temperature T are such that e~ 2/3 <C D <C 1. The first regime is dominated by the annihilation process 



G 



and leads to a configuration with less than three particles on the same site. This first regime occurs on a timescale 
of order 1. Then, the diffusion process comes into play on a timescale of order D~ x and the dimers diffuse until they 
reach a single particle and annihilate. Then the system reaches a configuration with mainly isolated particles. Finally, 
in order to reach the equilibrium concentration of particles, the activated regime involving the effective motion of 
isolated particles through the creation or annihilation of dimers is necessary and occurs on a timescale of order e 2 @ . 

The last regime in the concentration decay (before the equilibrium concentration is reached) may also be seen as 
either A + A^$otA + B^$ reaction-diffusion processes, depending on the models, since the particles have to 
pair themselves in order to disappear. Those two processes have different critical behaviours and critical dimensions: 
d c = 2 for the former and d c = 4 for the later. As a consequence, we expect a power law decay: 

c(t) ~ (e^/t) a (23) 

with a — 1 above the critical dimension d c , while below the critical dimension d < d Cl a — d/2 for the A + A — > 
case, and d/4 for the A + B — > case [14,15] (at the critical dimension there may be logarithmic corrections). 

In Fig. 4 (left) we present numerical simulations for different dimensions and a diffusive constant D = 1CP 4 of the 
concentration decay for the model with a single type of particles. The temperature after the quench is T = 1/10, so 
the different timescales are well separated (1 <C D^ 1 <C e 2/3 ), and the decay presents a two plateau structure. The 
first plateau is roughly independent of the dimension whereas the second plateau decreases with increasing dimensions 
to reach the mean-field value. Notice that the qualitative behaviour is maintained even in the mean-field limit. The 
dynamics during the last stage of the decay corresponds to a power law decay with the expected critical exponents 
a. Fig. 4 (right) shows similar results for the model with two types of particles, but with a diffusive constant D = 1. 
The diffusive timescale is now equivalent to the annihilation one, and the decay presents a single plateau structure. 
Again, we see the critical behaviour during the last stage of the decay, with the critical exponents a expected from 
the theory. 



(b) 




Time - Time 

FIG. 4. (Left) Concentration of particles for the model with a single type of particles after a quench from T = oo to 
T — 1/10 with a diffusive constant D = 1CP 4 and for different dimensions (d — 1,2 and 4 and mean-field). The A + A — > 
reaction-diffusion process during the second decay is illustrated by the change in power law for d < d c — 2. The two straight 
lines correspond to the expected exponents a — 1 and a = 1/2. (Right) Model with two different types of particles, and diffusive 
constant D = 1. The straight lines correspond to the expected power law decays corresponding to A + B — > reaction-diffusion 
processes in different dimensions: for d < d c — 4, a = d/4, while a = 1 for d > d c . 

In what follows we will concentrate mainly on the model with one kind of particle. 

IV. MEAN-FIELD SOLUTION 

The mean-field version of the model where all sites are neighbors allows to write exact evolution equations for the 
probability p l n (t) of a site i to be occupied by n particles at time t. The dynamical equations for p l n {t) read: 
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= -Pl>P3 -P3) - Dpi(p 2 + Ps ) + Dpi(p + Pl ) - e-^pl ) (l-p )+e-^p\p a , (24a) 

"rr = -PiP3 + P 1 P3 - Dp\{p 2 + P3 ) + Dpl(p +pi)- e~ 2 ViPo + e~ 20 plp o , (24b) 
dv i 

= -Pips + P\P3 - Dpi (po + Pl ) + Dpi (pa + p 3 ) - e" 2/3 p|p + e- 2 ^po , (24c) 

^ - -p 3 (l - p 3 ) +J4P3 - I>pj(Po + Pi) + £>pj(P2 +Ps) - e" 2/3 p 3 Po + e- 2 ^(l - Po). (24d) 

The time dependence has been omitted for conciseness and the sum over the neighbors has been performed using the 
fact that the probabilities pP n are indeed independent of the site j: 



N 

Notice that the right hand size of each of the equations (24) comprises three pairs of terms, each pair corresponding 
to a particular process, annihilation, diffusion and creation, with their respective rates. 

The dynamical equations satisfy the conservation of probability, ^2 n p l n = 1, which reduces the number of indepen- 
dent variables to only three, for example po,Pi and p2, in the mean-field model. The trivial equilibrium probabilities: 

P e n q = e- 0n (i:e-A (26) 
\fe=o / 

are a stationary solution of Eqs. (24). Equations (24) can be solved numerically to arbitrary accuracy. However, it is 
also helpful to consider their approximate solution, regime by regime, to better illustrate the underlying physics. 

A. Concentration decay after a quench 

After a quench from infinite temperature to a temperature T = 1/(3, the energy per site or concentration of particles 
c(t) decays from c(0) = 3/2 [p„(0) = 1/4 for all n < 3] to its equilibrium value 

c{^) = cM={j^ke-A(j^e-A . (27) 

\fe=0 / \fe=0 / 

If T is low and the diffusive constant small, such that e _2/3 <D<1, the energy decay may be decomposed into three 
different regimes. The first regime corresponds to the disappearance of sites with 3 particles due to the annihilation 
process and leads to a first plateau in the decay. The diffusion comes into play in the second regime and leads to a 
second plateau. On this plateau the particles are essentially isolated. The last regime is an activated one where the 
creation process is necessary to reach the equilibrium concentration. In previous models [11-13], where the diffusive 
timescale of dimers was equivalent to their annihilation timescale, the first plateau discussed here was absent. 

1. First regime: zero temperature and no diffusion 

In the first regime, we assume that only the annihilation process can occur. The equations (24) simplify, and only 
the two first terms on the right hand side are relevant. The probability p% decays to p% — 0, and, assuming that p 2 
stays constant, we get a good approximation for the different probabilities p n (t) (See appendix A and Fig. 5). The 
probabilities p n corresponding to the first plateau of the concentration decay are given by: 

Po- 0.458, pi~ 0.287, p 2 ~ 1/4, p 3 ~ 0. (28) 
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FIG. 5. Probabilities p n {t) during the first regime of the energy decay. From top to bottom: p ,pi,p2 and p$. Symbols 
correspond to simulations and lines to the numerical solution of the exact equations (full lines) and the analytical approximation 
(dashed lines). The difference between the two sets of lines is only visible for p2- 



2. Second regime: zero temperature and slow diffusion 



For slow diffusion, D<C 1, the first regime occurs on a timescale t ~ 0(1). For times t ~ 0{D~ 1 ) diffusion comes 
into play. The probability pz for a site to be occupied by 3 particles saturates to a quantity of order D. From (24d), 
noticing that dpz/dt ~ 0(D 2 ) and thus negligible, 

P^DJ^. (29) 
Po +Pi 

From the final probabilities p n from the first regime, given by Eq.(28), it is possible to determine the value of the 
second plateau in the concentration decay. Replacing the expression (29) for p 3 in Eqs. (24b) and (24c), we obtain: 

P- = (30) 
dp 2 1 - Pi 

which can be integrated from p~2 = 1/4 to P2(t) [with t ~ 0(D~ 1 )]. This leads to 

Inpi(t) -pi(t) - Inpi +pi = 2p 2 (t) - 2p 2 . (31) 
Knowing that the particles are essentially isolated during this second plateau, the probabilities p n (oo) = p n are then: 

Po-0.85, ft" 0.15, J5 2 ~0, p 3 ~0. (32) 



3. Third regime: activation and slow diffusion 

The third regime corresponds to the activated regime where the system has to overcome an energy barrier to reach 
equilibrium, and occurs on a timescale t ~ e 2/3 . An analogy in this third regime may be drawn with a reaction- 
diffusion process of the type A + A — > 0. When two isolated particles moving with an effective diffusive constant 
e~ 2/3 through successive creation-annihilation processes pair up, this pair diffuse on a faster timescale (<~ D^ 1 ) before 
they disappear when it reaches another isolated particle. We therefore expect a power law decay of the concentration 
during the last regime. 

In order to confirm this observation we first need to determine the probabilities p 2 and pz up to the order e~ 2/3 , 
neglecting dp2/dt and dp^/dt which are of order e~ 4/3 . From (24c) and (24d) we deduce: 

p 3 ^e- 2/ V D P 2~e- 2f3 {Dp + Pl ). (33) 
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Replacing these expressions in (24b), 



dpi 
dt 



we obtain the mean-field equation for the reaction-diffusion process A + A 
e~ 2/3 . The solution of this equation is: 



(34) 



with an effective diffusion constant 



Pi(*) 



Pi 



1 + 2e- 2 f%(t - t) 



(35) 



with t an initial time onto the second plateau of the concentration decay (D 1 <C t <C e 2/3 ). t ~ 501? 1 is a good 
estimate for D = 10~ 3 and T = 1/10 (see Fig. 6). 



From Eq. (35) it seems that the concentration does not reach its equilibrium value c e 



^ but decays to zero 



asymptotically. In order to account for the approach to equilibrium we need to include in (34) terms of order e~ 4/3 
previously neglected, which become relevant when pi <~ e~^. Taking into account leading order terms in e~ 4/3 (which 
means neglecting terms of order pie~ 4 ^, which are always negligible), we obtain: 



dpi 
dt 



-2e- 2 ? (p\ - *- 2 ' 3 



which has for solution: 



pi(t) = e 



p 1 + e-^-(p 1 -e-^e- 4e ~ 3 ^ t -~ t y 



(36) 



(37) 



For a timescale t ~ 0{e 2 ") we recover the behaviour (35) but for a timescale t ~ 0(e 3 ^) the concentration reaches 
the equilibrium one as: 



c(t)~pi(t)~e- (l + 2e- 4e ""(*-*)) 



(38) 



leading to an equilibration time T eq ~ e 3 " larger than the relaxation time r ~ e 2 " deduced from equilibrium autocor- 
relation. 
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FIG. 6. Concentration decay as function of time after a quench to the temperature T = 1/10 and for a diffusive constant 
D = 1CP 4 . Symbols correspond to numerical simulations, and lines to the analytical results for the first and third regimes. 
Inset: T = 1/6. 



B. Out— of— equilibrium correlation and response 



We now concentrate on the behaviour of two-time quantities, like correlation and response functions, in the out- 
of-equilibrium regime. 
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1. Correlation functions 



From the two-time out-of-equilibrium autocorrelation functions 

C n ,n'{t,t w ) = (S ni (t),ndni(t w ),ri) 

with initial conditions 

it is possible to construct all relevant two point autocorrelations, and in particular 

C(t,t w ) = {n t (t)n t (t w )) =^2nn'C n , n ,(t,t w ). 



(39) 



(40) 



(41) 



The correlation functions C nj „< correspond to the probabilities of having n particles at time t on a given site when 
there were n' particles at time t w < t on this particular site. They satisfy the equations (24) replacing p l n (t) by 
C n ,n'(t,t w ) but keeping the probabilities p n (t). This leads to an explicitly time dependent linear system of equations 
for the autocorrelations. 

In the following we concentrate on waiting times t w 3> D , that is, after the second plateau in the concentration 
decay. In this case, sites have mainly at most one particle. The probabilities to have 2 or 3 particles on a given site 
are small and satisfy the equations (33). The only relevant correlation function among C„ ; „/ is C\ t \, and C(t,t w ) ~ 
Ci.i(t,t w ). The correlation functions with n or n' larger than one are at least of order e -2 ^ smaller, whereas those 
with n = or n' = do not affect C(t,t w ). A simple calculation (see appendix B) leads the following differential 
equation for Ci,i 



dd 



dt 



— (t, t W ) — — Cl,l(t, ty 



l + 2pi(t) 



D 



1 + D 



-2/3 . 



1 



D 



1 + D 



Pi{t w )pi{t)e- 2f> . 



with the solution 



Ci,i(t, *„,) = pi(t) \ P i{t w ) + (1 - pi(t t0 ))e-( t -*™)/ T 



The correlation time r c is given by: 



1 + 



D 



1 + D 



(42) 



(43) 



(44) 



From this solution we recover that Ci t i(t, t) = pi(t), which is given by (37). 

Fig. 7 presents the autocorrelations Ci i) and C\^(t,t w ) for two different waiting times t w — 10 4 and 10 6 as 
function of the time difference t — t w . The temperature is T — 1/6 and the diffusive constant D = 10~ 2 . The analytical 
results Eqs.(37) and (43) agree with the numerical simulations for a system size N = 10 6 . 




,8 



FIG. 7. Out-of-equilibrium Ci,i(t,t w ) (circles) and Ci,i(t,t) (squares) after a quench to temperature T = 1/6. Symbols 
correspond to simulations and lines to the analytical result. Waiting times are t w = 10 4 (full lines) and 10 6 (dashed lines), and 
the diffusion constant D = 1CP 2 . 
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2. Response functions 



In order to determine the out-of-equilibrium response function, we need to introduce a perturbation at time t w after 
the quench. Different perturbations are possible but to get a response related to the autocorrelation, it is common 
to consider a random field on each site coupled to the corresponding observable [24]. The simplest possibility is to 
couple the random field to the single occupancy operator 6 nii i, leading to the perturbation 

SH = -h^dSn^. (45) 

i 

h is the strength of the field and will have to be small enough to stay in the linear regime. The random variables 
may be Gaussian variables or Ising spins (ej = ±1) with zero mean and unit variance. The corresponding (integrated) 
response function is the change in the expectation value of S n (t),i due to the perturbation, 

Xi{t,t w ) = h-^N-^Y,^n i {t),i)h (46) 

i 

where the overline stands for the average over the random field variables. This response is conjugate to the autocor- 
relation C\.i(t, t'), which as was shown above is the relevant one for long times and low temperatures, 

^E f ' £ i^wA(f),i) = C hi(t,t'). (47) 

The definition of the perturbation is not enough to determine the response. We also have to define how this 
perturbation affects the dynamical rules, maintaining the detailed balance conditions in order to ensure equilibrium 
asymptotically. Once again different definitions are possible. We will consider two of them. The natural definition of 
the perturbed dynamics is to use for the rates a Metropolis rule with the perturbed Hamiltonian H + SH, 

mm(l,e- l3A(H+SH ^ (M) (48) 

where A(H + SH) corresponds to the change in the perturbed Hamiltonian under the corresponding transition. One 
disadvantage is that this definition will only extract a response from unoccupied sites. A second possibility is to 
modify the dynamical rules by multiplying the unperturbed rates by another Metropolis factor 

min(l,e- /3A ^) x min(l,e- /3A ^) (MM) (49) 

It is easy to see that this modification of the dynamical rules preserves detailed balance with respect to H + SH. The 
advantage is that this definition allows to extract a response from occupied and unoccupied sites. For simple spin 
facilitated models the two dynamics yield equivalent responses, the second one being more efficient from the numerical 
point of view. We will see below that this equivalence does not hold for the present models. 

In order to determine the response function analytically we assume that only the site i is perturbed. Its probabilities 
Pn(t) are modified accordingly: 

Pn(*) =Pn(t) + he iX n(t,t w ) (50) 

where p n (t) are the unperturbed probabilities. The equations (24) have to be modified to take into account the change 
of rates. The zeroth order in h gives back the dynamical equations (24) for the unperturbed p n (t). The first order in 
h leads to a system of equations for the responses Xn(t,t w ) (see appendix C). 

For the perturbation (45), the only relevant response function is Xi (t,t w ). For the case of the modified Metropolis 
(MM) it satisfies (see details of the calculation in appendix C) 

!*?»>(«,<.) = - (i + *,w + ^) *rw + *.«> (i + jfj; - 0^) «-*. (.1) 

Neglecting the third term in the last parenthesis leads to the following approximation: 

X?°*(t,t») = (h>i{t)(l-e-W>) (52) 



12 



with the time t c already defined for the correlation function, Eq. (44). The timescale involved in the response and 
correlation function is thus identical. Fig. 8 shows the accuracy of this analytic result. The response function is 
non-monotonic, a common feature of the activated regime [9,25], and also observed in models of vibrated granular 
matter [28,29]. In the present case the non-monotonic behaviour is given by the fact that the response is the 
product of a decreasing function, pi(t), corresponding to the number of defects able to respond, and an increasing 
one, 1 — e ^( t ^*".)/ r = j corresponding to the monotonic rescaled equilibrium response function. The scaling from of 
the response function given in (52) is analogous to that in the Fredrickson- Andersen model [25]. The corresponding 
calculation for the case of the Metropolis (M) dynamics gives 

Xi M) (M™) - Xi MM) (M™) - 2/3A(t,t t0 ) (53) 

where 

A(t,t w )= Pl (t) [ Pl (t) -pi(^)e-(*-*»)/ r =] (54) 

The two responses (52) and (53) are different. Given that the second term inside the brackets in (54) decays faster 
than the first one [see Eqs. (35,37,38)], A(t, t w ) is always positive, and Xi^ < Xi MM \ as expected. 



(a) 



(b) 




t - 1,„ 



t - 1 , 



FIG. 8. Out-of-equilibrium response Xi(M™) f° r ( a ) MM dynamics, and (b) M dynamics, as a function of t — t w , at 
temperature T = 1/6, for waiting times t w = 10 4 (circles) and 10 6 (squares), and a diffusive constant D — 1CP 2 . The lines 
correspond to the analytical result. 



3. Fluctuation-dissipation relations 

Having obtained correlation and response functions, we can now study out-of-equilibrium fluctuation-dissipation 
(FD) relations [26]. Since we are considering the case of long but finite times, and therefore one-time quantities are still 
changing with time, FD relations have to be considered between the integrated response, Xi (t,t w ), and the difference 
of the conjugate connected correlation functions, c[ c \(t, t) — C[ c \ (t, t w ), where c[ c \(t, t') = C[ c \(t, t') — Pi(t)pi(t') (sec 
[25] for a discussion). In Figs. 9(a) and 9(b) we show the FD plots for the case of MM and M dynamics, respectively, 
for temperature T = 1/6 and waiting times t w = 10 4 and 10 6 (inset). 

There are several things to note. First, despite the fact that both response functions and difference of connected 
correlations are non-monotonic in t (which implies that the FD curves when plotted paramctrically for fixed t w start 
from the origin, go up, and then come back again), to a very good approximation Xi(i, t w ) — xi[C[ c l(t, t) — C[ c l(t, t w )], 
similarly to what was found for other simple strong glass formers [25] . Second, the FD curves approach the fluctuation- 
dissipation theorem (FDT) value as waiting time is increased, as expected. Third, the FD relations look almost linear 
(although this may be just a consequence of the fact that the departure from FDT is relatively small). In this case 
the FDT violation ratio X(t,t w ) [27,26] is just a function of the waiting time, X — X(t w ). 

We see from Figs. 9(a) and 9(b) that X > 1 for the case of MM dynamics, while X < 1 for the case of M dynamics. 
This can be traced back to Eqs. (43,52,53), which lead to T\i(t, t w ) = C'[ c \(t,t) — C'[ c l(t,t w ) ± A(t,t w ), where the 
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upper (lower) sign corresponds to MM (M) dynamics, together with the fact that A(t,t w ) > for all times. An 
approximation to the value of X(t w ) can be obtained from the following argument. The slope with which the FD 
curves leave the origin corresponds to the time regime in which the exponential in the second term of Eq. (54) decreases 
much more rapidly than pi(t) [see Eq. (35)]. If we assume that pi(t) does not change at all in this initial period, we 
may approximate A(t, t w ) ~ p\{t w ) (l — e~( t ~ tm ^ Tc ) . This in turn gives for the FD ratio X{t w ) ~ [1 =p Pi(t w )]~ , with 
the upper (lower) sign corresponding to MM (M) dynamics. For the plots of Fig. 9(a) this approximation predicts 
X(t w ) ~ 1.18, 1.05 for t w = 10 4 , 10 6 , while a linear fit to the data gives X m (t w ) = 1.11, 1.04 [X(t w ) ~ 0.85,0.95 and 
Xs. t (t w ) = 0.89,0.96 for Fig. 9(b), respectively]. 

Finally, in Fig. 9(c) we compare the behaviour in the mean-field model with that at finite dimensions. For d = 1 FDT 
is obeyed, similar to what happens in the Fredrickson- Andersen model [25]. For d > d c = 2 however, the FD plots 
coincide with the mean-field ones. This indicates that the aging behaviour is controlled by the out-of-equilibrium 
critical point of the underlying diffusion-annihilation process, and that mean-field serves as a good approximation 
for the physically relevant dimensions d— 2,3. 



(a) (b) (c) 




C c (t,t)-C c (t,t w ) C c (t,t)-C c (t,t w ) C c (t,t)-C c (t,t w ) 

FIG. 9. (a) FD plot for MM dynamics at T = 1/6 and waiting time t w = 10 4 (inset: t w = 10 6 ). The symbols correspond to 
simulations, the full lines to the analytical result, and the dashed line to FDT. (b) Similar plot for M dynamics, (c) FD plots 
in various dimensions (for MM dynamics): d = 1 (circles), d = 2 (squares), and MF (full line); T = 1/6, t w = 10 4 . 

V. CONCLUSIONS 

In this paper we have introduced simple models inspired by covalent glasses and topological froths. The models 
display strong glass forming behaviour despite their trivial thermodynamic properties due to the presence of constraints 
to the dynamics which generate dynamical frustration. We have studied the connection with underlying diffusion 
annihilation processes, and have shown that the aging dynamics of these models is dominated by the critical out- 
of-equilibrium fixed point of the associated diffusion-annihilation theory. We formulated a mean-field version of the 
models, which keeps most of the features of the finite dimensional ones, and which allows for an extremely accurate 
analytical solution of the low temperature out-of-equilibrium aging dynamics in the low temperature regime where 
activated processes play a dominant role. 
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APPENDIX A: FIRST REGIME IN THE CONCENTRATION DECAY 

The first regime in the concentration decay may be obtained by considering the dynamics at zero temperature and 
without diffusion (D = 0). The dynamical equations reduce to: 
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dpo 

dt 
dpi 

dt 
dp 2 

dt 
dp 3 

dt 



P3(Pl+P2) (Al) 

P3(Po-Pi) (A2) 

P3,{Pl-P2) (A3) 

-PziPo+Pi) (A4) 



with the infinite temperature probabilities p n = 1/4 as initial conditions. Notice that the right hand side of the 
equations is proportional to p% and the derivative of p 3 is negative indicating that a first plateau will be reached when 
p 3 = with non-zero values for the probabilities po,pi and pi- 

In the absence of diffusion, p2 remains approximately constant at its initial value, P2(t) — 1/4. Notice that if the 
initial configuration (before the quench) does not correspond to one at infinite temperature, this approximation is no 
longer valid. Using p + pi = 1 — P2 — P3 = 3/4 — p 3 and (A4), we are able to deduce: 

3e- 3 */ 4 

= 2 + e -3t/4 - ( A5 ) 

The probabilities po and p\ are determined considering the linear combinations a±po +Pi with a± = |(1 ± \/5) which 
satisfy the following differential equation: 

^ (a±Po +P\) = ^ (a±Po +Pi +a±/ 4 ) • (A6) 



A solution is given by: 



a ± + l + a 2 ± V2 + e- 3 */ 4 

We deduce the following approximations: 



4(a ± p (*)+Pi(*)+4/4) _( 3 \ 1/a ^ (A?) 



/,\ 1 3 + a/5 . V5~i 3 — a/5 , v^+i , . , 

^) = -4 + ^a7T^-^a7T^ ^ 

/ \ 1 1 + a/5 , vs-i a/5-1 , vs+i 

Mt) = -I + ^7T A 2 + ^7T A 2 (A9) 



with A = 3/(2 + e- 3 */ 4 ). 



APPENDIX B: OUT-OF-EQUILIBRIUM CORRELATION FUNCTION 

Here we determine the out-of-equilibrium correlations C n ,n'- They satisfy the equations (24), with the replacement 
of p l n by C n , n '- Since all correlations with n > 2 or n' > 2 are at least of order e~ 2 @, and the ones with n = or 
n' = are irrelevant for the correlation C(t,t w ), we concentrate on the set of equations corresponding to n' = 1: 



dCii 



-C AP3 + C 3 ,i(l -Ps) - £>G),i(P2 + Pa) + DC 2 ,i(po +pi) - e- 2 "C ,i(l -Po) + e- 2/3 C M p , (Bl) 

-C 1AP3 + C ,ip 3 - DC 1A (p 2 +P3) + DC 3 a(po+Pi) - e~ 2/3 C M p + e- 2/3 C 2 ,iPo, (B2) 

= -C2.1P3 + C hlP3 - DC 2 .i(Po+Pi) + DC .i{p2+P3) - e- 2fj C 2 ,ip + e^Cs.iPo, (B3) 

= -C 3 ,i(l - P3 ) + C 2 ,iP3 - DC 3 ,i(Po + Pi) + DCi,i(p2 +p 3 ) - e- 2/3 C 3 ap + e^Ou(l - Po). (B4) 

Equations (B3) and (B4) allow to obtain C 2j i and C 3j i as functions of C\x and Co,i, neglecting their derivatives and 
using (33), 
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DC 2 ,i = e- 2fi (C 0l i(pi +D) + Ci,ipi) , (B5) 

(l + D)C 3 ,i =e- 20 (Ci,i(pi+£>) + C o ,iPi). (B6) 

Combining all these results and the fact that Co t i(t,t w ) + C\ t \{t,t w ) —Pi(t w ), we deduce an equation for C\ t \\ 

^(t,t w ) = -C\i(t,t w ) (l + 2 Pl (t) + e- 2 " + J^j Pl (t) Pl (t w )e- 2 ". (B7) 

C(t,t w ) = pi(t)pi(t w ) is a trivial solution of this equation. A general solution is then Ci t \(t, t w ) = C(t,t w ) + C(t,t w ) 
with: 

^(t, t w ) = -C(t, t w ) (l + 2 Pl (t) + e" 2 ". (B8) 



The solution of this linear differential equation is: 



/ / / / 

C(t, t w ) = C exp 



(-^- | f 2e-\(t>t'). (B9) 



with t c = (1 + D/(l + L>))- 1 e 2 ' 3 . Using (34), we deduce: 



C(t,t w ) = C^-exp(- t —^). (BIO) 



Pi(t w ) 

The parameter C 1 is obtained from the initial condition C\^{t w ,t w ) — pi{t w ) leading to the correlation: 

Ci,i(t, t w ) = Pi(t) + (1 - pi(t t0 ))e-( t -*-)/^) . (Bll) 

APPENDIX C: OUT-OF-EQUILIBRIUM RESPONSE FUNCTION 

In this appendix we determine the out-of-equilibrium response function corresponding to a perturbation (45) in 
the case of MM dynamics (49). This prescription leads to a different response for spins with e = 1 and —1. The 
response is determined assuming only spin i is perturbed at a time t w after the quench at t = 0. The corresponding 
probability of occupancies are slightly modified compared to the unperturbed case: 

p\(t)=p n (t) + ke i x%(t,t w ) (CI) 

which defines the two sets of response functions Xn *t«) f° r a perturbed spin with different random variable ej = ±1. 
The total response for a given n is simply the average between the two responses for different e,. The response 
functions are determined from the first order expansion in powers of the field h of the modified equations (24) . \2 and 
X3 have a higher order in e~ 2/3 than xo and Xi- From the conservation of the probability, it follows that xo = — Xi 
and we may take Xi as the only relevant response function. In order to determine xi we need the equations for p\ 
and p\ for a random variable e, = 1 



dt 



dt 

while for a = — 1: 



-e-0 h p\p 3 + pip 3 - De-^ h p\(p 2 + p 3 ) + Dpl( Po + Pl ) - e^^iX) + e-^pjpo, (C2) 
-J4(l - Pa) +J4P3 - DpI(Po+Pi) + De-P h p\{p 2 +pz) - er^p^ + e^p^l - p ), (C3) 



dA 
dt 



= -P\P3 + e~ 0h pips - Dp\(p 2 +p 3 ) + De-e h pl(p Q + Pl ) - e~^p\p Q + e^^o, (C4) 



= -Pl(l - P 3 )+P 2 Ps - De-e h pi(p + Pl ) + Dp{(p 2+P3 ) - e-^plpo + e^p^l - p ). (C5) 
Using (33) and neglecting the time derivative of x|, the resulting equations for \n are: 
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^ = e-^(3 Pl (l+ Pl + D) - e-^ X +{l + 2 Pl + D) + D x + 



(C6) 



e ~2/3 



Xs+ = TTd ^ d ~ 0Pl {pi + ' (C7) 

^ = e- 2 ?f3 Pl (po + -D) — e-^xr (1 + 2pi + D) + D X * , (C8) 
P -2/3n 

Xs = TT ^-(xr-/3pi). (C9) 
From these equations we deduce closed linear differential equations for the response functions: 

d 4 =e-^(l + ^)-e-Xf(l + 2 , 1 + T ^), (CIO) 

= e - 2/? /3pi f 1 - Pi + T^) r- 2 "xi (l + 2 Pl + -^-) . (Cll) 



«it ^ V 1 + DJ /VL \ 1 + 

For systems where all sites i are perturbed with uncorrelated e,, by self-averaging and linearity, the total response 
X(t,t w ) - Xi(t,t w ) = (xf +Xi)/ 2 satisfies: 



dx 
dt 



(t, t w ) = (1 + 2px(t) + ^) X (t, ^K 2/3 + PMt) (i + ^ - ^§y) e " 2 "- ( C12 ) 



Discarding the third term in the last parenthesis allows us to solve this differential equation. This term may be 
neglected due to the small value of P i(t) for t > t w 3> D^ 1 or considering only the case of small diffusion constant D. 
Replacing x{t,t w ) — (3pi(t)(x(t,t w ) + 1) leads to the following equation for \: 

§(t,t w ) = -^A (CIS) 
at t c 

with t c the timescale already introduced for the correlation function, Eq. (44), r c = e 2/3 [1 + D/(l + D)] 1 . A solution 
satisfying the initial condition x(twjt w ) — or x{t w ,t w ) = —1 is: 

X(t,t w ) = - e -(*-*«)^« (C14) 

leading to the response function: 

X (MM) (t,t W ) = fait) (l - e -(*-*»)/^) . (C15) 
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